Spatio-temporal structure of cell distribution in cortical 
Bone Multicellular Units: a mathematical model 
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Abstract 

Bone remodelling maintains the functionality of skeletal tissue by locally coordinating bone-resorbing cells 
(osteoclasts) and bone-forming cells (osteoblasts) in the form of Bone Multicellular Units (bmus). Understand- 
ing the emergence of such structured units out of the complex network of biochemical interactions between 
bone cells is essential to extend our fundamental knowledge of normal bone physiology and its disorders. 
To this end, we propose a spatio-temporal continuum model that integrates some of the most important in- 
teraction pathways currently known to exist between cells of the osteoblastic and osteoclastic lineage. This 
mathematical model allows us to test the significance and completeness of these pathways based on their abil- 
ity to reproduce the spatio-temporal dynamics of individual bmus. We show that under suitable conditions, the 
experimentally-observed structured cell distribution of cortical bmus is retrieved. The proposed model admits 
travelling-wave-like solutions for the cell densities with tightly organised profiles, corresponding to the pro- 
gression of a single remodelling bmu. The shapes of these spatial profiles within the travelling structure can 
be linked to the intrinsic parameters of the model such as differentiation and apoptosis rates for bone cells. 
In addition to the cell distribution, the spatial distribution of regulatory factors can also be calculated. This 
provides new insights on how different regulatory factors exert their action on bone cells leading to cellular 
spatial and temporal segregation, and functional coordination. 

Keywords: bone cell interactions, cortical BMU, spatio-temporal bone remodelling, rank-rankl-opg, mathe- 
matical modelling 



1 Introduction 

In human adults, between 5 and 30% of bone volume 
is replaced every year [1, 2] in a process referred to 
as remodelling. Bone replacement is accomplished 
by stand-alone groups of cells of the osteoclastic and 
osteoblastic lineage progressing through old bone 
over a period of several weeks. Such a group of cells 
is called a "Bone Multicellular Unit" (bmu) and can 
be viewed as the basic functional unit for bone re- 
modelling [3, 4, 5, 6]. Tetracycline-based histomor- 
phometry has considerably helped in the elucidation 
of the spatial organisation and kinetic properties of 
the different bone cells in cortical bmus [7, 8, 9], 
which clearly indicates a spatial segregation of bone 
cell types depending on cell maturity. At the front of 
a bmu, in a region called the resorption zone (see Fig- 
ure 1), active osteoclasts attach to the bone surface 
and dissolve bone by secreting a mixture of proteases 
that break down the collagenous matrix, and hydro- 
gen ions that reduce the pH and dissolve the min- 
erals into the micro-environment [10, 11]. Towards 
the back of the bmu, in the so-called formation zone, 



active osteoblasts refill the cavity by laying down a 
collagen-rich substance known as osteoid, which sub- 
sequently mineralises to form new bone over the fol- 
lowing month or so (see [7, 1, 2]). The region be- 
tween the resorption zone and the formation zone, 
referred to as the reversal zone, contains precursor 
cells of both lineages embedded in a loose connec- 
tive tissue stroma [7]. New precursor cells and nutri- 
ents are supplied to the bmu by a small capillary that 
grows at the same rate as the bmu progresses into the 
bone. The net effect of the passage of a bmu at a spe- 
cific location of bone is the local renewal of the bone 
matrix and the formation of a so-called "secondary 
osteon", which includes a new Haversian canal. 

The existence of such a functional remodelling unit 
(referred to by Frost as a "packet of turnover" [3]) 
suggests the presence of tight couplings between the 
various cell types composing bmus. It has been hy- 
pothesised several decades ago that some combina- 
tion of local and/or systemic signals structure this 
internal cellular organisation [3, 6]. In the mid 
1990s, the discovery of the rank-rankl-opg pathway 
explained many previous experimental observations. 
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Figure 1 - Schematic figure of the internal organisation of a cor- 
tical bmu. Osteoclasts resorb the bone matrix at the front while 
osteoblasts lay down osteoid towards the back to refill the cavity. 
The central capillary provides a supply of precursor cells, as well as 
various nutrients. A schematic plot of the number of cells present 
at each position x along the bmu is depicted below. 



This regulatory pathway can be expected to play a 
major role in bmu physiology Many other potential 
regulatory molecules have been found by experimen- 
tal biologists (including systemic hormones, nerve 
signals, vascular agents, growth factors, chemokines, 
etc; see [12, 13, 14]). However, it is yet to be proven 
that these local interactions are able to group sev- 
eral generations of osteoclasts and osteoblasts in the 
form of bmus that present a clear spatial and tempo- 
ral separation of these cellular activities. While the 
structure of bmus is well understood at a descriptive 
level [7, 1, 2], how this structure is linked to the fun- 
damental underlying cellular interaction mechanisms 
remains to be elucidated. The present work aims to 
address this question. 

In this paper, we extend our previous tempo- 
ral model of bone remodelling [14] into a one- 
dimensional spatio-temporal model. Using this 
model, we study how bone cells structure themselves 
into a cortical bmu under the action of intercellu- 
lar signalling. This model is based on fundamental 
material-balance equations expressed as partial dif- 
ferential equations (PDEs). Non-conservative pro- 
duction or elimination of biochemical components in 
these general continuity equations are prescribed in 
accordance with the known biochemistry currently 
believed to play the most important role in bone re- 



modelling. As such, the model explicitly includes 
transforming growth factor ji (tgf-/3), parathyroid 
hormone (pth) and the receptor-activator nuclear 
factor k/3 axis consisting of the receptor rank, the lig- 
and rankl and the soluble decoy receptor osteopro- 
tegerin (opg). These regulatory factors couple two 
cell types of the osteoclastic lineage (a third one is 
introduced in Section 4) and three cell types of the 
osteoblastic lineage. Other components of the cel- 
lular communication system, known and unknown, 
are introduced implicitly through various model pa- 
rameters and external model conditions. For exam- 
ple, the capillary assisting the progression of a cor- 
tical bmu is modelled as a localised supply of bone 
precursor cells around the capillary's (growing) tip. 
Under these assumptions, we find that the model ad- 
mits solutions for the cell distributions in the form of 
travelling waves that have profiles that match the ob- 
served internal spatial organisation of a cortical bmu. 

In recent years, several teams of researchers have 
elaborated mathematical and computational models 
of bone remodelling, generally monitoring the evo- 
lution of the bone cells over time via ordinary dif- 
ferential equations (ODEs) [15, 13, 14]. Recently, 
Ryser etal. have included a spatial dimension in the 
model [15], addressing the important question of in- 
teraction between locally-expressed rankl and solu- 
ble opg for a trabecular bmu [16, 17]. In their model, 
bmus are driven by a rankl field in the surrounding 
bone matrix. Other researchers have developed cel- 
lular automata simulations to model resorption and 
formation on a per site basis [18]. 

To our knowledge, no group has yet addressed the 
issue of internal structuring of cortical bmus. Our 
approach emphasises the detailed integration of the 
biochemical processes involving osteoclastic and os- 
teoblastic cells at several maturation stages into a 
comprehensive partial differential model of the cor- 
tical bmu. Since it is based on a general formula- 
tion of the material-balance equation, the construc- 
tion of the model is modular and extensible. New 
interaction pathways or cell types can be included as 
needed. The one-dimensional continuous formula- 
tion employed here enables us to investigate analyt- 
ically how the various cell distributions making the 
internal structure of the bmu depend on the model 
assumptions. 

The paper is organised as follows: the model for- 
mulation is described in Section 2. In Section 3, the 
system of coupled nonlinear PDEs is then solved nu- 
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merically for the various cell and regulatory factor 
distribution profiles along the bmu. Theoretical inves- 
tigations of these profiles are performed, allowing us 
to map some of the profiles' properties to parameters 
of the model. In Section 4, we investigate the effects 
of various model assumptions made in Section 2. Fi- 
nally, we extend the initial model to include a new 
differentiation stage for osteoclasts, which is required 
to explain their observed spatial migration from the 
reversal zone to the resorption zone (see Figure 1). 
Concluding remarks are made in Section 5. 



since the fluxes are differential in space, they are ex- 
pected to play an important role in the spatial organ- 
isation of the cells within the bmu. 

In practice, the definition of local densities relies 
on a representative volume element large enough to 
contain many entities, yet small enough to remain 
local. While only few cells are present in a single 
bmu, continuous cellular densities can be defined in 
a statistical sense [20], i.e., by averaging histograms 
of cell counts over an ensemble of similar bmus (see 
Figure 1). 



2 Mathematical model of cortical 
bmu remodelling 

In the confined environment of a cortical bmu, the 
most important phenomena taking place are the bio- 
chemical interactions between the cells and their reg- 
ulatory factors, as well as the directed or diffusive 
transport of these entities. These phenomena are de- 
scribed in general by the material-balance equations 
of the species considered [19, 20, 21]: 

d 

—n A (r,t) = a A (r,t)-V-J A (r,t). (1) 
at 

In Eq. (1), A denotes any cell type or regulatory agent 
(such as hormones, growth factors, paracrine factors, 
etc.) explicitly accounted for in the model; n A (r, t) is 
the local density or concentration of A (i.e., number 
of entities A per unit volume) 1 at point r in space 
and at time t (r is a position vector); a A (r,t) de- 
notes local source/sink terms that account for non- 
conservative mechanisms, such as cellular prolifer- 
ation, differentiation, apoptosis, or mass action ki- 
netic rates of the regulatory factor binding reactions; 
J A (r,t) is the flux associated with transport prop- 
erties of A, such as diffusion, advection, or result- 
ing from inherent motility, e.g. chemotaxis. Due to 
the interactions between cells and regulatory fac- 
tors, the material-balance equations (1) written for 
all As are coupled. These couplings may arise both 
through the source/sink terms (e.g. hormonal up- 
regulation/down-regulation of a cellular response) 
and through the fluxes (e.g. chemotaxis). Note that 



: To align with common practice, we shall use the terminology 
"density" for cells and "concentration" for regulatory factors even 
if the units are chosen the same. 



Osteoclasts 

Following the ODE model of bone remodelling pro- 
posed by Pivonka etal. [14], we consider two stages 
of osteoclast development: "precursor osteoclasts" 
(oc p s) and "active osteoclasts" (oc a s). Precur- 
sor osteoclasts are assumed to have derived from 
hematopoietic progenitor cells and to be delivered 
to the bmu cavity at the tip of the capillary (see 
Figure 1). In cortical bmus, it takes 3 to 4 days 
for (single-nucleated) pre-osteoclasts to differentiate, 
migrate and join the dozen or so active multinucle- 
ated osteoclasts (each composed of around 10 nu- 
clei) found at the front of the bmu. These individual 
nuclei in active osteoclasts are then degraded after 
around 12 days [9, 7, 22]. In the model, oc a s repre- 
sent single nucleated entities incorporated in a mult- 
inucleated active osteoclast, and oc p s turn into oc a s 
upon RANKL-mediated activation of their rank recep- 
tor [11, 12]. Transforming growth factor /3 is known 
to be a general inhibitor of osteoclast differentiation 
and activation [11]. For simplicity, here we only as- 
sume that oc a apoptosis is enhanced by the presence 
of tgf-/3. Osteoclast maturation in the model can be 
summed up schematically as: 



We translate this sequence of events into the follow- 
ing balance equation for oc a s: 

d 

g^ n oc a = ®oc p ( RANKL )n c p - ^oc a ( TGF -rt "oc a ~ V - J oc a > 

(3) 

where @ oc is the RANKL-dependent differentiation 
rate of oc p s and j?/ 0Ca the TGF-/3-dependent apopto- 
sis rate of oc a s. As in Ref. [14], the up-regulation 
and down-reglation of cellular responses by a ligand 
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are assumed proportional to the fraction of occupied 
receptors. Mass action kinetics of the binding re- 
actions shows that this is equivalent to modulating 
the cellular responses by certain "activator" and "re- 
pressor" functions of the ligand concentration (see 
Refs. [13, 14, 23] for more details). With the dimen- 
sionless activator and repressor functions 

(4) 

the functional forms of @ oc and j?/ 0Ca can thus be 
written as 

% c (rankl) = D oc 7T act ( ^ ) , (5) 

oc p v J oc p \ RANKL I ' v J 

V^OCp J 

^(™-P)=A 0C y«(j^y (6) 

where fc"™ L and fc™ F "' 3 are dissociation binding con- 

OCp OC a 

stants, and D oc and A OCa are the maximal possible 
rates taken by @ oc and j4 ^ ■ 

In the confined space of a cortical bmu, cell diffu- 
sion is limited and we assume that directed motility 
dominates the movement of osteoclasts. The flux of 
oc a s can thus be written J OCa = n c a v oc a > where v OCa 
is the velocity of oc a cells with respect to the (fixed) 
bone matrix. The actual velocity of an active osteo- 
clast is a combination of the dissolution process of 
the bone matrix, and of chemotactic and/ or mechan- 
otactic signals [10, 11, 24, 25, 26]. Precisely how 
this sensing by osteoclasts of their mechanochemical 
micro-environment occurs is still uncertain and not 
an issue for the purposes of this paper. For this rea- 
son, in our model, the rate of movement of oc a s is 
simply taken to be constant, matching the average 
velocity u of the bmu's progression through bone: 

Joc a ="oc a "- (7) 

Note that typical cortical bmu velocities range from 
20 to 40 u.m/day [7, l, 2]. 

Osteoblasts 

Following the ODE model of bone remodelling pro- 
posed by Pivonka etal. [14], three stages of os- 
teoblast maturation are considered. "Uncommitted 
progenitor osteoblasts" (ob u s) denote a pool of mes- 
enchymal stem cells assumed to be provided around 



the tip of the capillary [8, 7, 1]. These cells are ca- 
pable of committing to the osteoblastic lineage, be- 
coming "pre-osteoblasts" (oB p s). This commitment is 
up-regulated by tgf-/3 [27, 28, 29] . Pre-osteoblasts 
further mature into "active osteoblasts" (ob 3 s), found 
in large numbers (1000-2000 cells) at the back of 
cortical bmus (see Figure 1), actively laying down os- 
teoid to refill the cavity opened by the osteoclasts 
[7]. Based on Pivonka etal. [14], osteoblast activa- 
tion is assumed to be down-regulated by tgf-/3. The 
fate of active osteoblasts is either to be buried in os- 
teoid and become osteocytes (approximately 95% of 
all bone cells are osteocytes), to undergo apoptosis, 
or to become so-called bone-lining cells covering the 
surface of the new Haversian canal [1]. Elimination 
of OB a s from the active pool is assumed here to in- 
clude all three possibilities. Osteoblast development 
in the model can thus be depicted as the sequence 

TGF-/3+ TGF-/3 — 

OB u > OB p > OB a > (8) 

leading to the following balance equations for OB p s 
and ob 3 s: 

d 

— t n 0B p = ®OB u (TGF-/3)n oBu - ®0B p (TGF-/3)n 0Bp - V-J 0Bp 

(9) 

d 

Yt n ° R - = ®ob p (tgf-/3) n 0Bp - A OBa n 0Ba - V-J 0Ba , (10) 

where @ 0Bu , @ 0B and A 0Ca are the ob u differentiation 
rate, the ob p differentation rate and the OB a elimina- 
tion rate, respectively. Similarly to Eqs. (5)-(6), we 
set 

® 0Bp ( TGF - rt = D 0Bp ^^j, (12) 
with k TGF '^ , k TGF ' 13 denoting dissociation binding con- 

OB u OBp 

stants and £> 0Bu , D 0B corresponding to the maximal 
possible rates taken by @ 0Bu and @ 0B . 

Active osteoblasts lay down osteoid in cortical 
bmus mainly radially, from the circumference of the 
cavity towards the center [8, 7, 1] . As this process oc- 
curs on a time scale much larger than that of resorp- 
tion, ob 3 s remain essentially stationary with respect 
to bone along the bmu axis. Furthermore, it is ob- 
served that active osteoblasts, unlike osteoclasts, are 
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not dynamically replenished once they have started 
producing osteoid [7]. This suggests that the pre- 
osteoblasts they derive from are not moving longitu- 
dinally either (at least, not to a significant extent), 
and so we set v OB = v OB =0, leading to 



0. 



(13) 



As will be seen in Section 4, this hypothesis is cru- 
cial to explain the spatial segregation of active os- 
teoblasts, pre-osteoblasts and uncommitted progeni- 
tors. 

Regulatory factors and binding reactions 

System-level coupling between the osteoclasts and 
osteoblasts occurs because the two direct regulatory 
factors in our model (tgf-/3 and rankl; see Eqs. (2) 
and (8)) are themselves driven by the cellular ac- 
tions, both directly and indirectly via other interfer- 
ing molecules. 

Tgf-/3 is stored in high concentration in the bone 
matrix and released into the bmu environment in ac- 
tive form by the resorbing osteoclasts [11, 28, 29]. 
Assuming that tgf-/3 degrades at a constant rate 
D TGF _fl, we have: 



dt 



n bone 

"TGF-jS ft -res' 1 oc a 



- V J T 



(14) 



where fc res is the bone volume resorbed per unit time 
by a single osteoclast and n'j™^ is the concentration 
of tgf-/3 present in the bone matrix. Since tgf-/3 
is released in the environment in soluble form, its 
transport properties encoded in J XGF _^ are assumed 
to be independent of the cells. It is expected that 
high levels of tgf-/3 are found up until the reversal 
zone where it promotes commitment and differentia- 
tion of mesenchymal cells to the osteoblastic lineage. 
For simplicity, we assume that tgf-/3 has negligible 
diffusion, i.e., J TG¥ _p «s 0. Nevertheless, the presence 
of tgf-/3 in the reversal zone can be accounted for by 
assuming a weak degradation rate D XGF _^ (in a sense 
clarified below). Further comments on the effects of 
tgf-/3 diffusion towards the back of the bmu are made 
in Section 4. 

The local availability of rankl, which is critical for 
the differentiation of oc p s into oc a s, arises from the 
combination of several effects. Rankl is a protein 
bound to the membrane of cells of the osteoblastic 



lineage. Its interaction with the rank receptor found 
on oc p is regulated by the presence of the soluble 
decoy receptor opg, which is also expressed by os- 
teoblastic cells. Furthermore, the relative expression 
of rankl vs. opg by osteoblasts is regulated by sys- 
temic pth concentrations. All these molecules and 
their competitive interactions are considered explic- 
itly in our model. Here we only describe their main 
features, and refer the reader to Ref. [14] for further 
details. We will assume that rankl is only expressed 
by OB p s and that opg is only expressed by OB a s (cor- 
responding to "Model Structure 2" of Ref. [14]). 
This choice of ligand expression is in agreement with 
experimental findings [30, 31] and the conclusions 
drawn in Ref. [14]. However, to reexamine this as- 
sumption in a spatio-temporal framework, we will 
study its influence in Section 4. While the flux of 
soluble opg is assumed independent of the cells (sim- 
ilarly to tgf-/3), transport of membrane-bound rankl 
is tied to the cells expressing it: J RASKL = n M(1 v 01 . 
However, osteoblasts are assumed to have negligible 
motility (v OB as 0), and so J RANKL ^ 0. 

A considerable simplification of the mass action ki- 
netic equations considered for the competitive bind- 
ings between rank, rankl and opg was obtained in 
Ref. [14] due to the separation of time scales be- 
tween the fast reaction rates of ligands binding to 
their receptors on cells, and comparatively slow cell 
responses. We examine here the consequence of this 
separation of time scales in the presence of transport 
terms in Eq. (1). Let r L be the slowest reaction rate 
(e.g., in day -1 ) to be found in the source/sink terms 
in a L for the ligand L. Dividing Eq. (1) by r L , one 
has 



dt 



(15) 



If reaction binding dominates transport, then 



\r7 V-J r | < 1 and rrV, 



0(1). Thus, changes in 



the local concentration of the free ligand occur on the 
short timescale r" 1 and only quasi-steady states need 
to be considered for the cellular dynamics, leading to 



i 



Vr,t. 



(16) 



This simplification is exactly of the same form as in 
the temporal model [14, Eqs. (16)-(20)]. We assume 
here that it holds for rankl, opg and for pth. As in 
Ref. [14], Eq. (16) is thus used to express the con- 
centrations of these regulatory factors in terms of the 
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remaining unknowns of the system. This has been 
presented in detail in Ref. [14], so we only briefly 
mention the results here. The pth endogeneous pro- 
duction rate P PTH (r, t) and degradation rate D PTH are 
assumed to be given and not further regulated. Thus, 
Eq. (16) with cr PTH (r, f) = P PTH — D PTH n PTH determines 
the pth concentration as 

^PTH -^PTh/^PTH (17) 

(see Eq. (25) of [14]). Production and elimination 
rates of rankl and opg in Ref. [14] have a more com- 
plicated form owing to their regulation by pth, the 
interdependence between rank, rankl and opg, and 
an assumed saturation of the endogeneous produc- 
tion responses. With similar notations as in Ref. [14, 
Eqs. (30)-(36)], the concentrations of opg and rankl 
can be rewritten with the help of the functions (4) as: 

„ _„ nr , „act f ft.°"c n °Bp + fe,QPG"0B a "\ 

"opg - OPG max? I ^— ^ 71 , 

V uru max 1/ opc J 

(18) 

_ AlANKL rep /, , \ 

^ RANKL p. *^ ^ Al, RANKL ^OPG * ^A2, RANKL ^RANK J 

RANKL 

„aCt f ^ RANKL ('t-jRANKL^ i r) RANKL „ ~\ _PTH \ r -\ (\\ 

X7Z o n on p +R 2 "oB a J7Iact,OB » (19) 

\ /"'rankl J 

In Ref. [14] and in the present model, the same con- 
stant number of rank receptors AT™ is assumed to 

OCp 

be expressed on each oc p cell. However, while the 
density of oc p s was constant in Ref. [14], it is space 
and time dependent here. The constant rank con- 
centration occurring in the Eq. (36) of Ref. [14] has 
thus to be replaced in Eq. (19) above by the local, 
time-dependent concentration 

n =JV RANK r! (20) 

"RANK OCp OCp • \.^>) 

Unlike Ref. [14], we do not assume that Eq. (16) 
holds for tgf-/3, and keep its differential description 
given by Eq. (14). Indeed, the production rate of 
tgf-/3 occurs on a cellular time scale and its degrada- 
tion rate is assumed to match the slower characteris- 
tic times of the cellular dynamics. 

We finally note that Eqs. (15)-(16) also apply 
to the balance of bound receptor-ligand complexes. 
Their fast binding properties allow us to express via 
(16) the receptor occupancy per cell in terms of the 
concentration of free ligand as has been used in 
Eqs. (5), (6), (11), (12) with the functions (4). 



External conditions 

Because all cells eventually differentiate further or 
undergo apoptosis, a continual supply of precursor 
cells is needed to reach nonzero cell populations over 
a period of time exceeding a couple of days. In corti- 
cal remodelling, this supply is local: it reaches the 
reversal zone of the bmu through an internal cap- 
illary that grows at the same rate as the bmu pro- 
gresses (see Figure 1) [7]. We assume here that the 
replenishment of ob u and oc p cells occurs around the 
tip of the capillary in an unbounded and non-rate- 
limiting way. Under that assumption, the inhomoge- 
neous densities n oc and n 0Bu instantaneously reach a 
stationary distribution peaked around the tip of the 
growing capillary [7]. These densities become given 
external functions in Eqs. (3) and (9), of the form 

"ocpCr, t) = oc p O-ut), rc 0Bu (r, r) = OB u (r-ut). 

(21) 

We assume oc p (r) and OB u (r) to be Gaussian dis- 
tributions centered around the capillary tip (see Fig- 
ure 2). 

Finally, while pth has been included into the model 
following Ref. [14], its spatial implications in the bmu 
will not be investigated for the purpose of the present 
study, and we assume that the concentration of pth 
is constant and homogeneously distributed along the 

BMU. 

Solving the system of PDEs (3), (9), (10) and (14) 
requires appropriate initial and boundary conditions. 
In the following, these equations are solved in one 
spatial dimension with boundary conditions specified 
at the very front of the bmu and at its back. 

3 Density profiles in the bmu 

As spatial profiles in a bmu are predominantly struc- 
tured along the longitudinal x-axis (see Figure 1), 
we restrict the mathematical model to this single 
spatial dimension, neglecting variations in transverse 
cross-sections: n A {r,t) » n A {x,t). As explained in 
Section 2, the fast binding approximation (16) al- 
lows to substitute rc PTH , n^, n RANKL , and n 0PG with 
their expression (17)-(19) in the PDEs (3), (9), 
(10) and (14), which are then solved numerically 
(using Mathematica [32]) for the remaining un- 
known concentration profiles n 0Ca , n 0B , n 0Ba and 
n TGF _p . These PDEs are of the reaction-advection type 
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and a boundary condition needs to be specified at a 
single point of the x axis in each equation. Based on 
bone physiology, we prescribe both n 0B and n 0B to 
be zero at the tip of the bmu cavity, and rjoth n 0Ca and 
n TGr .« to be zero at the back of the bmu, where the 
new osteon cavity is refilled with osteoid up to the 
diameter of the Haversian canal. These requirements 
in turn specify a bmu spatial domain over which the 
PDEs are solved. This domain is set on either side 
of the capillary tip (which moves along x at rate u) 
as follows. The tip of the bmu cavity is defined to be 
350 u.m ahead of the capillary tip while the back of 
the bmu is defined to be 4800 u.m behind the capillary 
tip [1], thus allowing the bmu to spread over about 
5 mm. 2 ' 3 To transform the moving-boundary condi- 
tions into time-independent conditions, the problem 
is solved in a reference frame co-moving with the bmu 
at rate u along x (see Ref. [33] for more details). 

3.1 Numerical results and discussion 

The evolution of the computed cell profiles is shown 
from the (static) bone frame in a series of tem- 
poral snapshots in Figure 2. These profiles define 
the shape of a multi-cellular wave front emerging 
and propagating into the bone at constant velocity 
u = 40 u.m/day, corresponding to a remodelling bmu. 
Starting at t = days from a tiny population of active 
osteoclasts present around the precursor cells (as- 
sumed to be recruited there during an "activation" 
phase of the bmu), the densities of ob p and oc a cells 
quickly increase to reach quasi-steady profiles at the 
front of the bmu over the next 20 days, progressing 
forward without changing shape. The tails of the ob p 
and (particularly) ob 2 profiles further at the back, 
however, take longer to develop. As a result of the 
differentiation sequences (2) and (8) the heights of 
the profiles automatically reach bounded values after 
an establishment phase and they are confined over a 
spatial range corresponding to the known dimensions 
of a bmu (of the order of a few millimetres). The de- 
velopment of a clearly structured travelling wave pro- 
file is predicted by the model, as is observed histolog- 
ically [7, 1, 2]. Pre-osteoblasts and active osteoblasts 
are shifted towards the back of the bmu. The inver- 



Evolution of the cell profiles (bone frame) 



Note, however, that cell densities can be concentrated on a 
more restricted portion of this domain. 

3 The origin of the (static) x-axis is also chosen such that it co- 
incides with the tip of the moving cavity at time t = 100 days (see 
Figure 2). 
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Figure 2 - Evolution of the cell profiles computed from the mathe- 
matical model. At t = days, the initial conditions are shown: 
only precursor cells and a tiny population of oc a s are present. 
At t = 3 days, these initial oc a s have released enough tgf-/3 in 
the environment to trigger differentiation of ob u s into OB p s, which 
in turn have increased the oc a population by RANKL-binding. At 
t = 20 days, the profiles at the front of the bmu have already taken 
a constant shape. Active osteoblasts have started to emerge be- 
hind OB p s. At t = 100 days, all profiles have essentially reached a 
quasi-steady-state: they are moving forward into the bone matrix 
without changing shape. A sketch of a typical bmu cavity is aligned 
with these steady-state profiles for comparison. Inset: OB p and OB a 
profiles represented in logarithmic scale together with the asymp- 
totic expressions given in Eqs. (28) (a = 0.036, b = -0.0545). 
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sion of the relative number of OB p s vs OB a s at around 
—850 [im (at t = 100 days) in Figure 2 shows that 
the model captures the transition between the rever- 
sal zone and the formation zone along the longitudi- 
nal axis of the bmu (compare also with Figure 1). 

On the other hand, it is apparent that the model 
does not capture the transition between the resorp- 
tion zone and the reversal zone: the populations of 
oc a s and oc p s overlap everywhere at the front of 
the bmu. The bone remodelling biochemistry impli- 
cated in the model thus far, whilst reproducing bone- 
formation features of the bmu very well, is not sat- 
isfactory in explaining the spatial cell structure in 
the resorption zone, which suggests there are miss- 
ing biochemical components not taken into account 
in this first model. In Section 4, we therefore supple- 
ment our model with a further cellular component to 
resolve this behaviour. 



3.2 Theoretical analysis of the cell pro- 
files 

The simple wave-form Ansatz n A {x, t) = A(x — ut) 
for the density profiles reduces the system of PDEs to 
a system of ordinary differential-algebraic equations 
(DAEs) to solve for the steady-state profiles A(x), A = 
OBp, OB a , oc a , tgf-/3 (see also Ref. [33]). 4 For OB p and 
OB a , these equations are 



- U ^-OB = % (TGF-0)OB U - @ (tgf-/3)0B 

OX p 

d 

- U ^-OB a = © (tgf-/3)0B - A OB a . (22) 

OX P Fa 

We can use Eqs. (22) to quantify the shifts of the 
osteoblastic profiles towards the back of the bmu as 
well as their relative height in terms of model param- 
eters linked to biological characteristics of the sys- 
tem. Let x* and x* be the positions of the max- 

OB p OB a f 

imum in the ob p and OB a steady-state profiles (in 
Figure 2, these are located at x* —550 jxm and 

OBp 

x* & —1000 \xm). Since the spatial derivative of 
ob p (x) vanishes at x* B , and that of ob 3 (x) at x* m , 
one obtains from (22) the following ratios of the den- 
sities of ob p vs. ob u and the densities of OB a vs. ob p at 



these points: 



-/(tgfP(x* b )), 



OB„ D 
— -fx* ) = — 



OB, ^ D 0B 

^«B a )=^Kw**: Ba )), 

UD p OB a 



(23) 



where 



V "-OB u y V "-OBp S 

g ( TG ^)=n^(l0) (24) 

are monotonously increasing and decreasing func- 
tions of tgf-/3, respectively. Due to the couplings 
existing between the various regulatory factors and 
the cells in the model, the tgf-/3 concentration oc- 
curring in the right hand side of Eqs. (23) depends 
implicitly on all cell densities (and in particular on 
their ratios ob p /ob u , OB a /oB p ), so Eqs. (23) do not 
express the cell density ratios in the left hand sides 
in closed form, and other dependences upon the 
fractions D 0B jD 0Bv and £ 0Bp /A 0Ba exist in / and g, 
respectively. Nevertheless, it can be checked nu- 
merically that the parameter fractions D 0B u /D 0B and 
D 0B /A m are main regulators of the cell density ra- 
tios in the left hand side of Eqs. (23). In fact, it can 
be argued that the implicit dependence of / and g on 
these parameter fractions enhances the explicit linear 
dependences seen in Eqs. (23). Indeed, upon increas- 
ing D 0B /D 0B , rankl is increased over opg, leading 
to an increase of oc a , thus of tgf-/3 and of /. On 
the other hand, upon increasing D 0B /A 0Ba , opg is in- 
creased over rankl, leading to a decrease of oc a and 
of tgf-/3, thus to an increase of g. 

Multiplying the second equation in (22) by 1 or 
by x and integrating it over the length of the steady- 
state bmu (from — oo to 0), one can use integration by 
parts, the boundary condition OB a (0) = and the fact 
that OB a (x) decays exponentially as x — » — oo (see Eq. 
(28)) to derive the following relations: 



r° 



dx % ob p -A OBs 



dx OB a , 



(25) 



r° 



dx OB a 



dxx@ 0B OB p -A OBa 



dx X OB,. 



4 The equation for oc a (x) becomes algebraical while the other 
equations keep a differential nature. 



Under the assumption that the variation of @ 0B along 
x can be neglected for the values of the integrals, 
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one can factor @ OB out of the integrals. 5 Given that 
the average position of the profile A(x) (or "center of 
mass" along x) is (x A ) = J dxxA[x)/ jdxA[x), one 
thereby obtains from Eqs. (25): 



Cell profiles with wrong cellular fluxes 



(*OB„) * (*OB„} - 



(26) 



Eq. (26) now shows that the ratio u/A 0B is essentially 
determining the length of the shift of the OB a profile 
towards the back of the bmu compared to OB p s. Per- 
forming a similar analysis with the first equation in 
(22), one has 



(27) 



meaning that OB p s are shifted to the back of ob u s by 
a length proportional to u/D 0B . 

Finally, it is possible to give analytical expressions 
for the decays of the ob p and OB a profiles at the far 
back of the bmu. Indeed, in this region, tgf-/3 is es- 
sentially absent, so that © 0B « 0, and @ 0B « D OB 
(see Eqs. (11), (12)). Eqs. (22) thus become linear 
and their solution can be calculated explicitly, leading 
to the asymptotic behaviours 



ob p (x) ~ a e |x| , 



(28) 



0Ba (x) ~ b e-^W + — ^- \e-^ - e-^Wl 



as x — » —oo, where a and b are two integration con- 
stants. These asymptotic behaviours of the steady- 
state profiles are compared with the numerical pro- 
files obtained by the temporal evolution in Figure 2 
in logarithmic scale. While the slope of these curves 
is essentially determined by the ratios -D 0Bp /u and 
A 0Ba /u, the constants a and b shift the height of the 
curves in the logarithmic plot, and they have been fit- 
ted. The small discrepancy visible at the very back of 
the bmu is due to the fact that the numerical profiles 
have not yet reached a quasi-steady state there. 

Equations (23), (26), (27) and (28) all relate ki- 
netic properties of the cells (velocity, differentiation 
and apoptosis rates) to intrinsic features of spatial 



5 In fact, @ OB and ® 0Bu vary significantly over the domain of 
integration. However, it is possible to use a generalised integral 
mean-value theorem (see [34, §27]) to factor these functions out 
of the integrals. This leads to exact relations between (x 0Ba ) and 
(x 0B > and between (x 0B ) and (x OB ) similar to Eqs. (26)-(27). 
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Figure 3 - Steady-state density profiles obtained when all cell ve- 
locities are taken equal, i.e., v 0Bp = v OBa = v 0Ca = u. All parameters 
are otherwise taken as in Figure 2. The typical spatial organisation 
of the bone cells in a bmu is not reproduced in this case. These pro- 
files are in clear mismatch with the typical shape of the bmu cavity 
sketched below. 



profiles in a cortical bmu at a fixed snapshot in time. 
This entanglement of time and space reflects the 
wave-like character of the bmu's progression. Impor- 
tantly, it is noted that the experimental observation of 
such profiles from histological analyses could allow a 
direct estimation of the cellular kinetic properties in 
this model. 



4 Roles of model assumptions 
and extension of the model 

In this section, we use our mathematical model to ex- 
amine the influence and assess the validity of several 
assumptions made in Section 2. Based on the discus- 
sion in Section 3.1 and the identified model short- 
coming, we then extend the model by adding a new 
stage of osteoclast differentiation. This extension re- 
solves the inability of the previous model to capture 
the transition between the resorption zone and the 
reversal zone. 

Influence of cell motility 

While the effects of many regulatory factors on cell 
commitment, proliferation, differentiation and apop- 
tosis are well-established in the context of bone re- 
modelling, less is known on the motile properties 
of the bone cells and the regulation thereof, al- 
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though recent progress has been made in this direc- 
tion [29, 25, 26]. 

Here we show that these motile properties can in- 
fluence dramatically the spatio-temporal coordina- 
tion of the bone cells, and thus the functioning of 
bone remodelling. In Section 2, OB p s and OB a s were 
assumed to stay stationary with respect to bone and 
we set their velocities to zero (see Eqs. (13)). The 
wave-like propagation of the density of ob p and OB a 
cells at rate u observed in Figure 2 is entirely due 
to their creation upstream and elimination down- 
stream. On the other hand, choosing ob p and OB a 
cell velocity to correspond to u, so that v 0B = v 0Ba = 
v 0Ca = u, still leads to a wave-like propagation of the 
cell densities at rate u. However, as can be seen from 
Figure 3, in this situation all the cell density profiles 
fall within the same region around the precursor cell 
source. Such cellular profiles are in clear mismatch 
with the experimentally-observed propagation of a 
structured cortical bmu with a separation in time and 
in space of the resorption and formation processes. 
Since all cells and regulatory factors overlap, their 
net interaction is modified and the cell density pro- 
files also reach different heights, offsetting the bone 
balance. This result corroborates the experimental 
observation that osteoblasts do not move significantly 
after their commitment to the osteoblastic lineage 
[7] and emphasises the critical role that cell motil- 
ity plays in BMU-based remodelling. 

Role of osteoblastic maturation stage for expres- 
sion of rankl and OPG 

Several experiments have characterised rankl and 
opg expression levels on maturing osteoblasts, con- 
cluding that rankl expression dominates in imma- 
ture osteoblasts while opg expression dominates in 
more mature osteoblasts [30, 31]. In the purely tem- 
poral model of Pivonka etal. [14], various "model 
structures" of expression of rankl and opg by os- 
teoblasts at different stages of maturation were 
tested, which supported these experimental findings. 
However, these model structures really need to be 
tested for their functional importance in the context 
of the spatio-temporal coordination of bone cells in 
a bmu. The density profiles predicted by our model 
when rankl is only expressed on ob 3 s and opg is only 
expressed on OB p s (corresponding to "Model Struc- 
ture 1" of Ref. [14]) are shown in Figure 4. Clearly, 
the lack of availability of rankl in the reversal zone, 



Cell profiles with rankl and opg expression swapped 
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Figure 4 - Steady-state density profiles obtained when the ex- 
pression of rankl and opg on precursor and mature osteoblasts 
is swapped, i.e., rankl is expressed on ob 3 s and opg is expressed 
on OBp. All parameters are otherwise taken as in Figure 2. The 
density of active osteoclasts has not grown past its small initial 
condition. Such a small population of oc a s would not be able to 
create a bmu cavity as before. 

due to its expression on OB a s at the back of the bmu 
and its further screening by opg expressed on OB p s, 
blunts activation of osteoclasts. Barely any oc a s are 
found in the steady-state. Such a tiny population 
of oc a s would not be able to create a bmu cavity 
with a size comparable to experimentally-observed 
bmus. These results thus strongly support the model 
structure where rankl is expressed on OB p s and opg 
is expressed on ob 3 s (i.e., "Model Structure 2" of 
Ref. [14]), which is used throughout the paper. 

Role of tgf-/3 

Tgf-/3 is a multi-facetted regulatory factor serving 
several purposes in bone remodelling [28, 29]. In 
Ref. [14] and in our model, tgf-/3 up-regulates os- 
teoblast commitment but down-regulates osteoblast 
activation (see Eq. (8)). Furthermore, it up-regulates 
active osteoclast apoptosis. These several roles find 
sense in the structured cell profiles of a cortical bmu. 
Physiologically, a negative feedback loop on osteo- 
clasts is needed to keep resorption under control. 
Having tgf-/3 regulating this negative feedback is 
convenient, since it is then activated only once bone 
has started to be resorbed. On the other hand, the 
portion of bone just resorbed needs to be refilled with 
new bone. While tgf-/3 diffuses behind oc a s in the 
cortical bmu, it reaches the reversal zone populated 
with ob u s. Tgf-/3 commands their commitment to 
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the refilling task by up-regulating their development 
into OB p s. Finally, the down-regulation of activation 
of ob p into ob 3 by tgf-/3 helps preventing osteoid to 
be laid down too early, e.g. in the resorption zone. 
From our model, we observe that the presence of tgf- 
/3 behind oc a s acts to delay the onset of OB a s and so 
to increase the spatial segregation between oc a s and 

OB^S. 



Cell and factor profiles — extended model 



RANKL + 

oc p — > oc r 



bone surface 



Or 



(29) 



To model the transition from oc m to oc a at the 
bone surface of the genuinely three-dimensional bmu 
cavity in our one-dimensional setup, we introduce a 
"probability of reaching bone" <£(x, t), defined as 

R\x, f ) 

#(x,t) = l-— ^-A (30) 



0J 

U 




o 

u 



Model extension: including "mature osteoclasts" 

Precursor osteoclasts are known to be circulating 

cells [12, 11] and delivered to the reversal zone of £ 

o 

cortical bmus through the capillary (see Figure 1). 
Throughout the bmu's progression, the capillary tip £ 
is found at a distance of about 350 Lim behind the re- 
sorbing front. This means osteoclasts need to travel 
this distance at a faster pace than the rate of progres- 
sion u of the bmu to reach the front [7] . On the other 
hand, activation of osteoclasts requires rankl, which 
is expressed on the surface of pre-osteoblasts found 
around the capillary tip. In the model presented in 
Section 2, oc a s are assumed to resorb the bone ma- 
trix, so while they have been activated by rankl, they 
cannot move faster than u. For this reason, oc a s are 
differentiating from oc p s around the middle of the 
reversal zone and stay there as they have not been 
given regulatory mechanisms to distance themselves 
from their progenitors (such as chemotactic signals 
towards the bone surface [26]). This results in over- 
lapping oc p and oc a density profiles in Figure 2. 

To resolve this problem, we are led to introduce 
a new stage of osteoclast development in the model, 
that we call "mature osteoclasts" (oc m s). Mature os- 
teoclasts denote osteoclasts that have been activated 
by rankl and that migrate towards the front at a 
velocity v 0Cm > u until they reach the bone surface. 
Once at the bone surface, they join an active resorb- 
ing multinucleated osteoclast, hence becoming oc a s 
progressing at rate u. The sequence of osteoclast mat- 
uration in Eq. (2) is thereby extended to 
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Figure 5 - Cell density and regulatory concentration profiles in a 
cortical bmu as given by the extended model (compare with Figure 
2). Active osteoclasts are now clearly shifted towards the front of 
the bmu compared to their progenitors. The overlap between OB a 
and oc a is practically inexistant. The osteoblastic profiles in the 
back are essentially unchanged except for their amplitude. The 
concentrations of opg, rankl and tgf-/3 are also shown. They have 
been scaled to correspond to ob 3 , OB p and oc a density levels, re- 
spectively for comparison. Remarkably OB p -bound free rankl is 
not colinear with the presence of OB p : it is effectively shifted to- 
wards osteoclast precursor cells due to the presence of opg. The 
concentration of tgf-/3 is also found behind the oc a peak. 



where R(x, t) is the radius of the bmu cavity depicted 
in Figures 2-5, and R c = 100 rim is the maximal cav- 
ity radius (cement line radius). 5 For simplicity, this 
cavity is assumed here to be given and to progress 
forward at rate u without changing shape (an as- 
sumption valid in the quasi-steady-state). According 
to Eq. (30), the function 3>(x, t) thus interpolates be- 
tween one ahead of the (moving) cavity front (where 
the cavity radius 7?(x, t) is formally zero), and zero 
in the reversal zone (where R(x, t) reaches the ce- 
ment line radius R c ). The function t) increases 
again to a positive value < 1 for x in the formation 
zone (where R{x, t) decreases towards the Haversian 
canal radius). With this definition, the sequence of 



This probability corresponds to the bone volume fraction of a 
cross-sectional slice at x of a rotationally-symmetric cortical bmu 
of radius R(x, t) compared to the cylinder of radius R c . 
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osteoclast development in Eq. (29) can be expressed 
as: 



'oc™ ^oc„ 



®oc>ankl) n oc - -% c ($)n 0Cm - V-J oc , 



dt 
d 

g^ n oc a = ®oc m O)"oc m ~ ^oc a (TGF-/3)n OCa - V-J 0Ca , 

(31) 

where @ oc , j4 oc?l and J 0Ca are given by Eqs. (5), (6) 
and (7). The transition from oc m to oc a is assumed 
to take place at a rate ® 0Cm proportional to the prob- 
ability of reaching bone: 



(32) 



where D nr is the rate of osteoclast activation once at 
the bone surface. Finally, the flux of oc m s is taken to 
be 



(33) 



where the velocity is such that v 0Cm > u. 

The simulation results of this extended model are 
presented in Figure 5. While the cell density profiles 
at the back of the bmu have not changed qualitatively, 
the front of the bmu now also exhibits structured pro- 
files, clearly delineating a transition from the resorp- 
tion zone (predominantly populated with oc a s and 
oc m s) to the reversal zone (predominantly populated 
with precursor cells). This structure reproduces the 
histologically expected cell profiles of a cortical bmu 
as schematically depicted in Figure 1. 

The concentration profiles of opg, rankl and tgf-/3 
are also plotted in Figure 5, with their heights scaled 
to correspond to the maximum of the ob 3 , ob p and 
oc a density profiles, respectively. One can see that 
tgf-/3 has slightly diffused to the back of oc a s and 
that its decline towards the back coincides with the 
onset of OB a s. As a result, there is no overlap of ob 3 s 
with oc a s. Since transport of soluble opg has been as- 
sumed slow compared to its reaction processes (see 
Eq. (16)), opg is found mainly in the same location as 
OB a s that express it. On the other hand, while rankl 
is bound to the membrane of OB p s, the rankl concen- 
tration profile is clearly shifted towards the front of 
the ob p density profile, overlapping in particular with 
oc p and oc m cells. This is due to opg inhibiting most 
of the available rankl ligands at the back of the ob p 
profile. With this observation, the role of opg takes 
on a new fundamental meaning. 



5 Concluding Remarks 

In this paper, we have developed a spatio-temporal 
mathematical model of cortical bmu remodelling 
based on fundamental material balance equations ex- 
pressed for different bone cell types. This model is an 
extension of the purely temporal model of Pivonka 
etal. [14]. Tgf-/3, the rank-rankl-opg pathway and 
pth are explicitly taken into account in the model, 
and mass action kinetics is used to describe the reac- 
tion rates between ligands and their receptors. The 
resulting system of (nonlinear) PDEs is solved for the 
cell densities and regulatory factor concentrations 
in one dimension, corresponding to the longitudinal 
axis of a cortical bmu (see Figure 1). 

We find that the cell interaction pathways in the 
model are able to explain the emergence of a multi- 
cellular travelling wave that develops structured pro- 
files moving without changing shape. The spatial 
structure of these steady profiles corresponds to the 
known organisation of bone cells in a typical corti- 
cal bmu. It clearly delineates a resorption zone at 
the front, followed by a reversal zone, and a forma- 
tion zone at the back. Several properties of the cell 
density profiles have been linked theoretically to ki- 
netic properties of the cells in the model, such as dif- 
ferentiation and elimination rates. These rates may 
thus be directly inferred from the measurement of 
cell counts in serial histological sections taken at a 
particular time point. 

It is experimentally known that several maturation 
stages of osteoclasts and osteoblasts have different 
expression patterns and responses in the tgf-/3 and 
rank-rankl-opg pathways. Our model shows that 
this heterogeneity is essential to ensure a functional 
BMU-remodelling process with segregated but coordi- 
nated zones of resorption and formation, in particu- 
lar: 

• Tgf-/3 plays a central role in modulating cell re- 
sponses as soon as bone is resorbed. It mod- 
erates osteoclastic resorption and initiates os- 
teoblastic commitment of mesenchymal cells 
once it has diffused from the resorption zone to 
the reversal zone. Furthermore, high concen- 
tration of tgf-/3 towards the front of the bmu 
helps prevent osteoblasts from becoming acti- 
vated prematurely, or near active osteoclasts. 

• The fact that rankl is bound to the membrane 
of pre-osteoblasts helps ensure that osteoclasts 
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do not initiate resorption without the presence 
of bone-refilling precursor cells. By shielding 
the availability of free rankl on OB p s, the ex- 
pression of soluble opg by active osteoblasts in 
the formation zone provides a mechanism to 
"shift" the peak of the concentration profile of 
free rankl towards the front of the bmu, ahead 
of the peak of the ob p population that expresses 
rankl, and thus prevents activation of osteo- 
clasts where new bone is being laid down. Our 
model shows that changing the rankl and opg 
expression pattern fails to coordinate oc a forma- 
tion properly. 

• The various motile properties of bone cells are 
absolutely crucial for the spatial organisation of 
the cells into a cortical bmu, both in the forma- 
tion zone and in the resorption zone. In partic- 
ular, our results suggest that osteoclasts migrate 
forward at various rates depending on their mat- 
uration, and corroborate the observation that 
osteoblasts are quasi-stationary with respect to 
bone. The importance of bone cell motile prop- 
erties is expected to play a critical role in the 
balance between bone resorption and bone for- 
mation, both in health and disease. 

While our one-dimensional model is capable of re- 
producing important features of cortical bmus, there 
are a variety of possible improvements that could 
be addressed using the general formulation of the 
model presented in this paper. Solving the system 
of PDEs in higher spatial dimensions could be used 
to study how cells distribute themselves in transverse 
cross-sections of the bmu. Other specific cell interac- 
tion pathways could be included as needed and stud- 
ied in detail on the basis of the present model. We 
note that investigating initiation and termination sig- 
nalling mechanisms of cortical bmus is very important 
to fully understand bone homeostasis during remod- 
elling, and will be the subject of future work. 
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